
#include <Eigen/Dense>
#include <fstream>
#include <iomanip>
#include <iostream>
using namespace std;

struct Fun1 {
  double coef{1};
  double operator()(const double &x) const { return x * coef; }
  double div(const double &x) const { return coef; }
};

struct Fun2 {
  double coef{1.0};
  double operator()(const double &x) const { return x * x * x * coef; }
  double div(const double &x) const { return 3 * x * x * coef; }
};

struct Fun3 {
  double coef{1.0};
  double operator()(const double &x) const { return x * x * x * x * x * coef; }
  double div(const double &x) const { return 5 * x * x * x * x * coef; }
};

// void init()

int main() {

  double x[4]{-0.2, -0.1, 0, 0.1};
  double DX = 0.1;
  int interp_order_ = 4;
  Eigen::MatrixXd vad(interp_order_, interp_order_);
  double ix = -0.01;
  for (int j = 0; j < interp_order_; ++j) {
    double dx = x[j] - ix;
    cout << "dx = " << dx << std::endl;
    vad(j, 0) = 1.0;
    for (int k = 1; k < interp_order_; ++k)
      vad(j, k) = vad(j, k - 1) * dx;
  }
  std::cout << "vad = \n" << vad << std::endl;
  Eigen::MatrixXd inv = vad.inverse();

//   vad /= DX;
//   for (int i = 0; i < interp_order_ - 1; ++i)
//     inv /= DX;

  double div{0};
  Fun2 f;
  for (int i = 0; i < interp_order_; ++i)
    div += f(x[i]) * inv(1, i);
  double exact = f.div(ix);
  cout << setprecision(16) << scientific;
  cout << "exact = " << exact << "  div = " << div
       << "   delta/exact = " << (div - exact) / exact << endl;
  return 0;
}